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We investigate quantum degenerate Bose gases at fi- 
nite temperature in the grand canonical ensemble [l[ and 
use the Bose-Hubbard model which can describe the dy- 
namics of ultracold atoms in periodic potentials such as 
optical @, i, H and magnetic lattices [1, S 0, II 1, E3 ■ 
A superfluid of ultracold bosons trapped in a periodic 
lattice undergoes a superfluid to Mott insulator quan- 
tum phase transition if the barrier height between the 
lattice sites is adiabatically increased @, S HI- I n the 
Mott insulator phase there is a fixed number of atoms 
per lattice site which may have applications in quantum 
computation [2j, |3j . The Bose-Hubbard model is also im- 
portant for the study of systems with strongly correlated 
bosons (Tlj . 

Ultracold bosons at zero temperature in the Bose- 
Hubbard model have been studied usin g di fferent meth- 
ods such as Monte Carlo simulations fl2l fl3l [T3 . [Hj], 
mean field theory QJ, Q3, GJ], Bethe-Ansatz solu- 
tion [l9j], exact diagonalization [2(J, strong-coupling 
expansions [2lj |. densit y- m atrix renormalization group 
(DMRG) (infinite-size) @, DMRG (finite-size) [23j and 
exact diagonalization plus renormalization group [2~j| . 
DMRG [23[ gives high precision results in only one di- 
mensional many-body problems [25l . |2(| . Bose-Hubbard 
model at finite temperature has been considered via 
mean- field theory [lfl [27], |28|. Monte Carlo simulations 
of a quantum rotor model [29(, an ab initio stochastic 
method [30l |. p erturbative DMRG [3l| and slave particle 
techniques [321 ] . So far, most of the finite temperature 
studies on the Bose-Hubbard model have been based on 
perturbation theory or some approximations. 

In this paper we use gauge P representation which is 
an exact phase space method based on a coherent state 
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representation. We promote the use of gauge P represen- 
tation as an exact method to benchmark various approxi- 
mate methods and simplifying assumptions. We evaluate 
the performance of this method for these imaginary time 
calculations of the Bose-Hubbard model to calculate cor- 
relations at finite temperature to connect between the dif- 
ferent limiting cases. Applying a phase space method for 
a system with a Hamiltonian written in the second quan- 
tized form, it is possible to convert the master equation 
of the system (such as the Liouville-von Neumann equa- 
tion) into a differential equation [33[ . Phase space meth- 
ods such as positive P representation [§3] and gauge P 
representation (35[ can give accurate results which their 
accuracy depends on the number of simulation trajecto- 
ries. Using the gauge P representation, the second-order 
spatial correlation function and also momentum distri- 
bution have been calculated for an interacting ID de- 
generate Bose gas in the Bose-Hubbard model We 
investigate ultracold atoms at finite temperatures with 
the gauge P representation [35[ and open boundary con- 
ditions. We choose the gauge P representation over other 
quantum phase space methods including the positive P 
representation [34J because in studying the many-body 
physics problems with strongly correlated bosons, the 
gauge P representation gives more stable results and also 
the sampling error is reduced, compared with the posi- 
tive P representation [3(1 Furthermore, the positive 
P/gauge P is exact, whereas other phase-space methods 
are not. 

We have verified the gauge P technique by compar- 
isons with exact numerical and also analytical results in 
simple cases [38[. Our simulation results are, within the 
sampling error, in remarkable agreement with the known 
limiting cases when either the hopping matrix element 
J or the onsite interaction U is zero. Therefore, they 
could be considered as a touchstone for testing the re- 
liability of approximate techniques when both J and U 
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are nonzero. Nevertheless, with the present gauge in the 
gauge P representation, because of increase in sampling 
error, simulation results for the the physical quantities in 
the Bose-Hubbard model are not precise at temperatures 
below T = 0.05U, when both J and U have finite values. 

In ID, we simulate ultracold atoms in the Bosc- 
Hubbard model with up to 11 lattice sites and study 
the average number of particles and coherence between 
lattice sites at finite temperatures. We show that for 11 
lattice sites, the edge effects are not very important and 
for a Bose-Hubbard model with 11 sites, we calculate the 
Luttinger liquid parameter which is important for locat- 
ing the boundaries of the Mott insulator lobes (39l. I40ll4l| . 

The Bose-Hubbard Hamiltonian is @, 0] 
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where M is the number of lattice sites and <i,j> 
means that the summation is taken over adjacent sites 
only. Also, aj, Sj and n% — a\a.i are creation, annihi- 
lation and number operators, respectively. The canon- 
ical commutation relations for fij and &t are [oj,ttJ] = 
5ij. The hopping matrix element J is defined by J = 
- / rf 3 x«i*(x - x 4 )[-?i 2 /2toV 2 -I- Vo(x)]w(x - Xj), where 
Vb is a periodic potential like the optical lattice (see [42j , 
for a review on optical lattices) or the magnetic lat- 
tice d, H, [13] potential. Wannier functions w(x — Xf) 
are localized position eigenstates 43] . The on-site inter- 
action U — g J d 3 x|ui(x)| 4 , where g — 4:Tra s ti 2 /m and a s 
and m are the s-wave scattering length and mass of the 
bosonic atom, respectively. 

At zero temperature, if U^>J, the system is well into 
the Mott insulator (MI) regime. For the Mott insulator 
with the commensurate filling of rij = n = N/M, where 
rii are the number of particles per lattice site, the ground 
state of the system for J = is described by | y , f = 
\n, ■ ■ ■ ,1%) (44|. Atom-atom correlations Ci(r) |22l. l23l |45| 
and density-density correlations Di (r) [45] are defined by 

Ci(r) — /a]d^_)_ r \ and Di(r) — {hihi +r ), where r is an 

integer and < r < M. The standar d deviation of the 
number of particles [g3| is Arii = \J (nf) — {hi) 2 . For 
the ground state given by \^ MI ) ]=a = \n,--- ,n), the 
atom-atom correlations, the density-density correlations 
and the standard deviations at each site are Ci(r) — 
0, Di(r) — n 2 and Arii = [45[, respectively. For a 
Mott insulator with J =^ 0, according t o first-order per- 
turbation theory, we have An.; = J /U \J 'Adn(n + 1) [46| . 
where d and n are the dimension of the system and the 
average number of particles in a Mott insulator lobe, re- 
spectively. 

At zero temperature, when the on-site interaction U 
is very small compared with the hopping matrix element 
J, the system is a superfluid (SF) and the ground state 
for U = can be given as a coherent state |^ SF ) [/=0 = 
y/ N\/M N £ {M |m, • • • , n M )/Vnil • • • n M \ [H, where 
J2i n i — N is the total number of particles. Also, the 



sum extends over all sets of occupation numbers {ni} 
subject to < ni < N. For M ^> 1 and commensurate 
filling N/M = 1, we obtain d{r) = 1, A(r) = 1 and 
Anj = 1 [45]. Therefore, for the commensurate filling 
of n = 1 we have the same density-density correlations 
Di(r) — 1 for both the superfluid and the Mott insulator 
ground states. 



I. FORMALISM OF THE GAUGE P 
REPRESENTATION AND DERIVATION OF iTO 
STOCHASTIC EQUATIONS 

In this section, we perform some calculations in the 
positive P representation and finally write the Ito form 
of the Langevin equations in the gauge P representation. 
Then, in section |H] we give the simulation results. Also, 
in this section to obtain some general expressions valid 
for both cases where U ^ and U = 0, using the Boltz- 
mann constant fcs = 1, we define dimensionless imagi- 
nary times t = U/T for U ^ and r' = J/T for U = 
and J ^ 

Now, we consider the unnormalized density matrix 
p u = e -(. Hr -t lNT ) _ which is in the grand canonical ensem- 
ble, and take its derivative with respect to r = U/T [l| 



dp u 
dr 
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where [^4, B] + = AB + BA is the anticommutator of A 
and B. Also, p e is defined by 



We now have 

dpu 
dr Z 

where 
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According to Eq. (|B1[) . in the positive P representation, 
we have 



P(a,/3,r)Ad 4M A 
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Taking the derivative of Eq. © with respect to r and 
considering Eq. (|4]) and Eq. (|6]), we obtain 



dP(a,P ,r) 
dr 



Ad iM X = 



1 J p(a,/3,T)H'(a, at)Ad 4M A 
-i J P(a,p,T)kH'{a,aJ)d 4M \ 
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We have, from Eq. (|S4| . H' N (a,aJ)A = H' A {a,j3 
d a )A and AH' N (a 7 a 1 ) = H' n (ol + dp,/3)A where 



M M M 
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Ad A = / P(a,/3,r)£^ j Ad 4 ^A (10) 



where 



Choosing the gauge gu = i\J%{n' k — \nk\) [35[ and con- 
sidering Eq. (fr7|). we can now calculate the Ito equa- 
tions and convert them into the Stratonovich form of 
the Langevin equations (see appendix [C] for more calcu- 
lations details). The Stratonovich stochastic equations, 
which are a natural physical choice |47j and more suitable 
for numerical solutions, compared to the Ito stochastic 
equations, and also have superior convergence proper- 
ties (sEEi, are 
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According to Eqs. ([8} and (|9]), we also have 
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Comparing Eq. (|13l) with Eq. (IA7|) . we obtain 
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where the Wiener increments dWi have the property [34l . 

m 

(dWi(T)dW^(«)) a = %<S(r - s)dr 2 (20) 
where () , means stochastic average. 



II. SIMULATIONS AND COMPARISONS IN 
DIFFERENT LIMITING CASES 

We simulate the ID Bose-Hubbard model Stratonovich 
equations Eq. (Tl9l) using extensible Multi-Dimensional 
Simulator (XMDS) [64[. We use XMDS with the semi- 
implicit interaction picture SUP method. We assume 
that at r = 0, there are N particles on average, in a 
thermal state [331 ] in the system. So the initial conditions 
are 
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(15) 

Dtj = -^%« j2 , t,j = 1,2,-- • , 2A7 (16) 

where w, 3 = tf-y-i + ^t-ij, ^i+Mj+M = ^ji and 
^>i+M,j = ^i,j+M — 0. Also, the effective complex boson 
number is rij = + = aift = n i+ M- Considering 
Eq. (|B8| . we have 



(17) 



* 3 - = vw/2(& + *e OT ), ft- = a*,, Q = 1 
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where ^ and ^ m are random numbers with zero mean and 
standard deviation 1 which can be realised by Gaussian 
random numbers with zero mean and standard deviation 
1 1331. 



A. One-site model 

Figure [T] compares the highly accurate numerical cal- 
culations, using a truncated number state basis, with 
the gauge P representation simulations for M = 1 [Hj]. 
It shows the independence of the expectation value of 
the number of particles (hi) at the target temperature 
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FIG. 1: Solid lines show the gauge P numerical calculations 
of (hi) versus r for M = 1, J = and no = 0.5, 1, 2,3 and 
4.5. The target chemical potential \i T is 0.5. Solid, upper 
dotted, and lower dotted lines show simulation results (hi), 
(hi) + a (sampling error) and (rii) — a, respectively. Black 
dashed lines show the highly accurate numerical calculations 
using a truncated number-state basis. 



T = 10 U (here, r T = 0.1) from the initial average num- 
ber of particles no, which is in good agreement with 
the numerical calculations based on a truncated number- 
state basis. As the figure shows, for a sufficiently large 
target imaginary time r T (sufficiently small temperature 
T) and a common value of the chemical potential \x T , (ni ) 
is independent of the initial number of particles no [37| . 
According to the figure, the expectation value of the num- 
ber of particles at r T = 0.1 is 2.66, which is independent 
of the different initial values of no = 0.5, 1, 2, 3 and 4.5. 
Independence from the no is very useful for the phase 
space simulations, because for different sets of the chem- 
ical potential and other parameters such as J and U the 
sampling error and also the stability of the simulations 
depends on no- 

Figure [5] shows the average number of particles, (ni), 
versus the inverse temperature r for the large value of 
t t =20. There is good agreement between the highly 
accurate numerical values and the gauge P simulation 
results. Figure [3] shows the simulation results for (hi) as 
a function of the chemical potential u T . By decreasing 
the temperature sampling error is increased, specially at 
high values of the chemical potential. It is possible to 
reduce the sampling error by increasing the number of 
simulation trajectories. According to Fig. [3J the stepwise 
pattern of the average number of particles as a function 
of the chemical potential vanishes as the temperature is 
increased from T = 0.1 U to T = 10 U. Because there is 
no hopping matrix element J, this result is valid for an 
array of M lattice sites in ID, 2D and 3D. Therefore, for 
J = 0, in a system with M lattice sites in ID, 2D and 3D, 
the stepwise pattern in the (hi)-fj, T plot vanishes as the 
temperature is increased from T = 0.1 U to T = 10 £/. 
For the Bose-Hubbard model, according to Eqs. (JSJ) and 




2.5 



& 2 



1.5 



(b) 




FIG. 2: (fit) versus r for M = 1, n = 3, J = 0, U = 1.0, 
dr = 10 -3 , u T = 0.9 (fi e = 0.9144) and (a) n p = 2 x 10 9 and 



(b) n p = 2 x 10 10 . Convent 
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B. Two-site model for U = 

For M = 2, we have a double well and the Stratonovich 
equations can be written from Eq. (|C13j) and Eq. (|19[) . 
Figure 2] shows (n.2) for a double site model (M = 2) 
with U = 0, J = 1, no = 5 and fi T /J = —1.01 
(He/J = —0.9918), where u r = —1.01 is the target chem- 
ical potential at T = 0.1 J. This figure also compares 
the exact analytical results [38J with the gauge P simu- 
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FIG. 3: Gauge P simulation of (ni) versus fi T , for A/ = 1, 
no = 1.2, J = 0, [/ = 1 and three different values of r T . 
Here again, fi T is the chemical potential at the target tem- 
perature T. As the temperatures is increased from T = 0.1 U 
to T = 10 U, the stepwise pattern of the average number of 
particles versus the chemical potential, which is clearly seen 
below To = 0.06 U as a truncated number-state basis shows, 
vanishes. Dots around the diagram for r T = 10 show the sam- 
pling error. For T — U and T — 10U sampling error is too 
small to be seen. 
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FIG. 4: Comparison of the simulation results for (712) with 
the analytical results, for no = 5, M = 2, U = 0, J = 1 
and n T I J — —1.01 {fJ. e /J = —0.9918). Here, we have r' = 
J/ksT — 1/T and the number of the stochastic trajectories 
n p is 10°. Other conventions as in Fig. [T] 



lations and shows good agreement between them. The 
sampling error is due to the stochastic initial conditions 
given by Eq. ([2lj). 



C. Two-site model for U / and J =£. 

The simulation codes have been tested for the limit- 
ing cases when either U = or J 7^ and are now 
ready for the general case of a two-site system in which 
both hopping matrix element J and the on-site interac- 
tion U exist. In Fig. [5] (h 2 ) (blue solid line), the rel- 
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FIG. 5: Simulation results for the expectation values of 
the number of particles (712), the relative standard devia- 
tion An2 = An2/ (712), the relative standard deviation for 
a coherent state with the same number of particles An co h = 
An co h/ {712) and the relative standard deviation for a thermal 
state with the same number of particles An t h = An t h/ (^2) 
at one of the sites in a double well system, r is proportional 
to the inverse temperature and the dotted lines around each 
line show the sampling errors. Here, the stochastic averages 
have been taken over n p — 10 9 trajectories and the target 
chemical potential fi T is 0.5. 



ative standard deviation A712 = An 2 / (^2) (black solid 
line) , the relative standard deviation for a coherent state 
with the same number of particles An co h = An co h/ (A2) 
(green dashed-dotted line) and the relative standard de- 
viation for a thermal state with the same number of par- 
ticles Anth — Anth/ {fi2) (red dashed line) in a dou- 
ble well system are shown, where An2 = y (n|) — (ni) 2 , 

An™/, = yj (h 2 ) and An th = J (h 2 ) 2 + (h 2 ). Here the 
on-site interaction U is 1 and J/U in Fig. 03(a) and (b) 
is 0.1 and 1.0, respectively. In this quantum simulation, 
the number of trajectories, n p , is 10 9 . Also, the target 
chemical potential, fx T , is 0.5. 

As Fig. El[a) shows, when J/U is 0.1, the average num- 
ber of particles at the temperature T = 0.05 U (r T = 20) 
is almost 1 which is close to the exact value of the aver- 
age number of particles for J/U — 0, shown in Fig. ??. 
Note that, as we increase J/U to 1.0, at r T = 20, we have 
(77,2) ~ 1.8. Also, if we increase the hopping matrix ele- 
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FIG. 6: (h) and the relative standard deviation An = An/(n) for the hopping matrix element J = 0.4 and three different 
system sizes M — 3, 7 and 11. Here, we also have no = 1.2, U — 1.0 and u T = 0.5. When the size of the system changes from 7 
to 11, for the central sites, 4 and 6, respectively, the expectation values of the number of particles and especially the standard 
deviations, within the sampling error, are in good agreement. 




ment the relative standard deviation An 2 increases and 
in both cases for J/U = 0.1 and J/U — 1.0 we have 



An 2 < An coh < An t h (23) 

It is interesting also that the relative standard deviation 
An 2 approaches that of a coherent state as we move to 
large values of J/U which at low temperatures are in 
the superfluid regime. This is in the right direction for 
describing a superfluid as a coherent state. 



III. SIMULATION OF M-SITE MODEL FOR 

17/0 AND J ^ 

For M > 3, the Stratonovich equations can be written 
from Eq. (|C13|) and Eq. (Il~9l) . Figure [6] compares (n) 
and An = An /(ft) at the central site for three different 
system sizes with the hopping matrix element J = 0.4. 
This figure shows that the size effect is not very consider- 
able for M — 11 as, within the sampling error, the values 
of (n) and An are similar to those for M = 7 and, in 
particular the relative number fluctuations are the same. 
Therefore, a Bose-Hubbard model with M = 11 may be 
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aia 2 ) I (fie) and ( asfig) = (fisSg) / (rie) at two different target temperatures T = U in (a) and (c) and T = 0.5 [7 in (b) 



and (d). The labels in (a)-(c) are the same as in (d). By decreasing the target temperature (increasing r T ) Coh, Ke, ( &1&2 



and (^a^al^ increase. Also, at a constant temperature, by increasing the hopping matrix element, all of these measures of the 
coherence between lattice sites increase. 



large enough to approximately represent the behaviour of 
an infinite system. In Fig. [7] the relative standard devia- 
tions Ahq, , are plotted versus the inverse temperature, 
t, for a general state and are compared with those of an 



ideal thermal state Anth, where Anth 



and a coherent state An co h, where yAn co h = 

two different values of J/U . In order to obtain An for 
the thermal and the coherent states we have used the 
simulation results for (h). According to the plots, at a 
constant low temperature, as J/U is increased, the stan- 
dard deviation increases. Also the standard deviation at 
the target temperature T and the target chemical poten- 
tial /x T for both different values of J/U is less than the 
standard deviation for a coherent state and much less 
than that of a thermal state. Also because at low values 
of J/U the atoms form a superfluid, when the tempera- 
ture goes to zero, the closeness of the standard deviation 
An@ to An co h supports the idea of describing a superfluid 
by a coherent state. Also an increase in the (quantum) 
standard deviation (n) is observed as the temperature 
increases. 

Because the number of atoms changes with tempera- 
ture, it is important to know the behaviour of the mea- 
sures of the coherence relative to the average number of 
atoms per lattice site. Figure [5] shows relative coherences 



«i4) = (ai4) / («e) and ^d 5 d^ = (asa^) I 



and 



also Coh = Coh/ (fie) and Ke = Ke/ (rig), where Coh 
and Ke are defined as 



Coh 



( ; r > a U 

\M(M- 1) ^ 1 



M-l 



Ke = \ 2(M-1) ? + 



(24) 



(25) 



Coh and Ke show coherence between all lattice sites and 
between all the adjacent sites, respectively. According to 
Eqs. (UJ and (|2"4"|) . between the expectation value of the 
kinetic energy part of the Bose-Hubbard Hamiltonian 



M-l 



KE = -J Ua\a i+1 + a l a\ +li 



(26) 



and Ke, Eq. (1251) . which is an average of the coherence 
between adjacent cites, there is a simple relation 



-KE 

Ke=^r^- r 7T, M^l 



2J(M-1)' ~ ' * (2?) 
As Fig. [5] shows, the relative coherence between sites 
I and 2, (dial ) , is different from that between sites 5 
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and 6, (a^Og), due to the edge (size) effect. Because the 

edge effect in a system with size M = 11 is small and 
all the other adjacent sites, except sites 1 and 2, have 
almost the same coherence (or same relative coherence) 
between each other as exists between sites 5 and 6, we 
have, according to Fig. [5J at temperatures T\ = U and 
T2 = 0.5 U (within the sampling error) 

Ke = (asa\) (28) 

According to Eq. ([2"5"j) and Fig. |5J for a Bose-Hubbard 
model with M=ll sites in ID, at temperatures 7\ = U 
and T 2 = 0.5 U, the kinetic energy of the system is simply 

KE = -20j/a 5 4)- 

In general, for M lattice sites, the kinetic energy term 
KE, Eq. (HD, is given by -2J(M - 1) (a M/2 a\ M+2)/2 ^ 

for even M and —2J{M— 1) ^a(M-i)/2^(M+i)/2^ ^ or oc ^ 
M > 3. 

In Fig. [8j all the parameters show an increase when 
either the system is cooled or the hopping matrix ele- 
ment J is increased. By decreasing the temperature at a 
constant J the relative coherence increases; therefore the 
ultracold bosons approach a superfluid state where the 
coherence between the lattice sites is very high. When 
the constant temperature is low enough, decreasing the 
tunnelling rate takes the system to a Mott insulator phase 



where coherence between the sites is lost [2j, |3( . Figure 
IHlshows (fie) and the relative values Auq, (^a^a^J, Coh, 

J(he), JAri6, J^a^a^j and JCoh as a function of J for 

two different temperatures. Here again, according to the 
figure, all the parameters increase when either the tem- 
perature is reduced or the hopping matrix element J is 
increased. 

Figure [10] shows C-zir) = C*2(r) / (ho) as a function of t. 
All values of C2(r) show an increase when T is reduced. 
The same behaviour is observed when J is increased. 

According to Figs. [8] and [TOl all relative measures of 

the coherence Coh, Ke, (aia^, (a^a^) and Ca(r) in- 
crease when either the temperature is reduced or the hop- 
ping matrix element is increased. 

IV. CALCULATION OF THE LUTTINGER 
PARAMETER 

The important parameter which has been used to lo- 
cate the critical parameters of the superfluid to Mott 
insulator quantum phase transition is atom-atom cor- 
relations C(r) = ^aja r ^ [H, [23[ and more generally 

Ci(r) = (a\a i+1 ^ against J [Hj]. 

According to Refs. [H, H(| 52, basically, it is possi- 
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FIG. 10: Relative atom-atom correlations Ci{r) = Cii?)! i^) as a function of r for M=ll, fi T = 0.5, values of J = 0.3 and 
0.5 and two different target temperatures T = U in (a) and (c) and T = 0.5 U in (b) and (d). The labels in (a)-(c) are the 
same as in (d). The relative atom-atom correlations increase when either the temperature is reduced or the hopping matrix 
element is increased. 



1.5 



: 



5~ 



M = 0.5 



0.5^ 
J 

0. L 



0.1 




k J = 0.3 



1 



1.5 



5~ 



0.5 



T 2 = O.hU 



K = 1.6 ± 0.1 



\J = 0.5 



| \ "* ♦ 

t V K = 3.8 ± 0.2 
J = 0.3 *i. 

H=P- 1 K = ^ i ±6:^ i> . 

;— ■■■■■■^ ...... 



4 5 

r 



4 5 
r 



(b) 



FIG. 11: Atom-atom correlations C*2(r) as a function of r, where r is the site number, for M=ll, /i T = 0.5 and two different 
target temperatures (a) T = U and (b) T — 0.5U and values of J = 0.1,0.3 and 0.5. Diamonds and circles around them 
show the simulation results and their sampling error, respectively. Dashed lines, which look like dotted lines, are fits to the 
function r~ K ^ 2 where K is the Luttinger parameter. According to this figure, C2(r) increases as temperature is decreased. 
Also at a constant temperature by increasing the hopping matrix element J, coherence between sites increases. Moreover, as 
J is increased the Luttinger parameter K decreases. 



ble to fit the function r~ K l 2 to the curve of Ca(^) as a 
function of r for different values of the chemical potential 
and hopping matrix element, find the boundaries of the 
Mott lobes, and obtain the phase diagram of the Bosc- 
Hubbard model. Figure [TT] shows the atom-atom corre- 
lations Ca(r) as a function of r for T — U and T — Q.5U. 
At each temperature three different values of the hop- 
ping matrix element J are considered. According to this 
figure, C'2(r) increases as the temperature is decreased. 
Also at a constant temperature by increasing the hopping 



matrix element J, the coherence between sites increases. 
Moreover, as J is increased the Luttinger parameter K 
decreases. At fi T /U = 0.5 and at the constant temper- 
atures T-y = U and T 2 = 0.5 U, by increasing the hop- 
ping matrix element J, relative atom-atom correlations 

C^(r), where C 2 {r) — {^'Wr^ which are other mea- 
sures of the relative coherence between sites, increase but 
the Luttinger parameter K, which is important in locat- 
ing the boundaries of the Mott-insulator lobes, decreases. 
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Also, for a constant value of the hopping matrix element 
,/ = 0.5, by reducing the temperature from T± = U to 
T2 = 0.5 U, the Luttinger parameter K changes from 
4.0 ±0.2 to 1.6 ±0.1. 



V. CONCLUSION 

In this paper, using the Bose- Hubbard model, we simu- 
lated ultracold atoms in the grand canonical ensemble of 
quantum degenerate gases. We studied ultracold atoms 
at finite temperatures with the gauge P representation. 
We have written a simulation code using XMDS, which 
generates a C++ code. In ID we simulated the Bosc- 
Hubbard model with 1, 2, 3, 7 and 11 sites. 

The simulation results are in good agreement with 
highly accurate numerical calculations, based on a trun- 
cated number-state basis, when there is no on-site inter- 
action between atoms. Also, for a double well system, 
we compared the simulation results with exact analytical 
results for the case where the atoms can tunnel between 
sites but the on-site interaction between the sites is zero. 

We also investigated the average number of particles, 
relative standard deviations and coherences between sites 
at finite temperatures for the Bose-Hubbard model in ID 
with open boundary conditions consisting of 1, 2, 3, 7 and 
11 sites and showed that at non-zero temperatures the 
relative standard deviation is not zero even for J/U = 
and grows as J/U is increased. We found that the relative 
standard deviation An^ is higher at higher temperatures 
and is less than the corresponding relative standard de- 
viation for a coherent state and is much less than that of 
a thermal state with the same value of (hi). 

For J = 0, we showed that above T — 0.1 U the step- 
wise pattern in the plot of (hi) versus the target chem- 
ical potential fj, T starts to vanish; therefore, there is no 
Mott insulator-like lobe in the phase diagram of the Bosc- 
Hubbard model in the plane of (i/U-J/U. 

Comparing the Bose-Hubbard model with M — 3, 7 
and 11, we showed that in a ID Bose-Hubbard model 
with M = 11 and open boundary conditions edge effects 
are insignificant except for the side sites (the first and 
last sites). 

At low temperatures, for constant values of J/U = 0.2 
and \i T /U = 0.5, by reducing the temperature from 
T\ = U to T2 = 0.5 U, the relative standard deviation 
of the number of particles at each lattice An^ decreases 
but remains well below those of a coherent state, An co hi 
and also a thermal state, Anth, with the same average 
number of particles. The relative standard deviation An; 
is closer to An co h than to Anth- This confirms that the 
ground state of a superfluid can be described by a coher- 
ent state. 

For n T /U = 0.5, at the constant temperatures T\ = U 
and T2 = 0.5 U, the relative standard deviation An^ de- 
creases as J/U is reduced. At these temperatures the 
lowest value of the relative standard deviation of the 
number of particles at each lattice Arii is not zero even 



if the tunnelling is zero. This confirms that at these 
temperatures, for J/U = 0, the compressibility K is not 
zero, so there is no Mott insulator phase present at these 
temperatures. This is in good agreement with the tem- 
perature T = 0.1 U discussed above or, more accurately, 
with the melting temperature To = 0.06 U above which 
there is no Mott insulator phase in the Bose-Hubbard 
model |38|. 

At the constant temperatures T\ = U and T2 — 0.5 U, 
for fi T /U = 0.5, by increasing the tunnelling rate all 

physical quantities Coh, Ke, /aiajV (0-50^ and 62(7), 

which measure relative coherence between the lattice 
sites, increase. Also, if J/U is kept constant but the 
temperature is reduced from T\ = U to T2 — 0.5 U, all 
of them increase. Likewise, all measures of coherence 
increase when either the temperature is reduced or the 
hopping matrix element is increased. 

For a Bose-Hubbard model with M=ll sites in ID, at 
temperatures T\ = U and T2 = 0.5 U, the kinetic en- 
ergy part of the system is simply given by —20 J ^ effigy, 
within the sampling error, which can be generalised 
to KE = —2J(M — 1) (aM/2 a \ M+2 )/2) ^ 0T an even 
number of lattice sites M, and to KE = — 2J(M — 

1) (cl(m— i)/2^(Af-|_i)/2 / f° r an °dd number of lattice sites 
M. 

At u T /U = 0.5 and at the constant temperatures 
T\ = U and T2 — 0.5 U, by increasing the hopping ma- 
trix element J, the Luttinger parameter K, which is im- 
portant in locating the boundaries of the Mott-insulator 
lobes, decreases. Also, by reducing the temperature from 
T\ = U to T2 — 0.5 U, for a constant value of the hop- 
ping matrix element J — 0.5 the Luttinger parameter K 
changes from 4.0 ± 0.2 to 1.6 ± 0.1. 

With the particular gauge-choice used here, we have 
found that the the growth of sampling error is a limiting 
factor at low temperatures and a large number of sites. 
However, it is always possible that for particular situa- 
tions a better choice of gauge may reduce the sampling 
error. 

Even with the sampling error limitations described 
above, the gauge-P method is very general and could 
potentially be applied to a range of ultracold lattice sys- 
tems, including 

1. The J/U critical ratios and phase diagrams in ID, 
2D, and 3D, for the Bose-Hubbard model and also 
the critical values and phase diagrams of the super- 
fluid to Mott insulator quantum phase transition at 
zero temperature. 

2. Disordered Bose-Hubbard model @, H3, [EH and 
also two component bosons in periodic lattices [5l[ 

3. The coexistence of the Mott insulator and su- 
perfluid phases in inhomogeneous traps, such as 
quadratic and quartic trapping potentials, for a 
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continuous range of incommensurate fillings [54|, 

4. Ultracold bosons in a double- well potential [5a| and 
a tilted multi-level double- well potential [57| . 

5. Extended Bose-Hubbard models [58| . 

6. Strongly interacting bosons in a 2D rotating square 
lattice which can also be studied via a modified 
Bose-Hubbard Hamiltonian I59I1. 



When there are no terms higher than second order, C A 
may be expanded as 



(+) 



(A7) 



In the positive P representation, the Fokker-Planck equa- 
tion is given by 

^gf' T) = (V-djA^ + ^djD^P^frr) (A8) 



Appendix A: Positive P representation 

There are density operators for which the P represen- 
tation does not exist. For example, the P representa- 
tion cannot describe nonclassical effects such as squeez- 
ing or antibunching [33[. In contrast, the positive P rep- 
resentation gives stochastic differential equations which 
can represent genuine quantum-mechanical problems like 
squeezing or antibunching (6^, HH, HH . 

In the positive P representation, we can write [34| 



where A 
and 



p = J P(a,(3,T)Ad 4M \ 
(a,/3) = (ai,a 2 , • • • ,a M ,0i,02, 
|a}</3* 



A = 



(/3*|a) 



= \\t X )(f3*\\e- a -e 



(Al) 

,0m) 

(A2) 



1 1 a) is an M-dimensional Bargmann coherent state [331 ] 
and \\ai) — e' Qi '/ 2 |cti), where M is the number of modes 
(number of lattice sites, for example), r may represent 
real time or imaginary time (inverse temperature) and 



exp 



(A3) 



where i — 1, 2, • • 
aA = 

Aa = 



, M. We also have the identities 



aA, a* A = 



^ da 



A, 



d_ 

d/3 



A, Aa f = /3A. 



(A4) 



Expectation values of the normally ordered products 
aj m a" can be written in the positive P representation 
as the following 



- fm " r 

ai a] 



(A5) 



Also, after integration by parts, provided the boundary 
terms vanish at infinity, we can write 



dP(a,/3,T) 



Ad 4M A= / P(a,/3,r)4 +) Ad 4M A 



(A6) 



Appendix B: Gauge P representation 

In the gauge P representation, the density matrix can 
be written as 35] 



where a 
(Q, a, (3), and 



p = J G{d,r)kd AM+2 c 



(a°,a\--- ,a M ,a M+1 ,a M + 2 , 



(Bl) 



„2M\ 



a = n 



\<*)((3* 
(P*\<x) 



= n\\a)(^*\\e- a ^ 



(B2) 



When fl = 1, this phase space representation reduces 
to the positive P representation j35|. Here, we have the 
identities given by Eq. (1A"4|) plus fi^A = A. 

In this phase space representation, quantum averages 
of the normally ordered products aj m a™ are 



at m a? 



stoch 



(n + n* 



(B3) 



stoch 



where (f (a)} stocl} = J f(a)G(a, r)d± M + 2 d. 

After integration by parts, provided the boundary 
terms vanish at infinity, we have 



dG( ^ T) Ad iM a= I C'{ n. t)C,. a Ad LW+ ~n (LU) 



In the gauge P representation, when there are no terms 
higher than second order, C GA may be expanded as 



C - C {+) 

GA '-'A 



V + 7^9 ■ ffidn + g k B jk dj 



(Vd n - 1) 
(B5) 



C 



OA -= + nDd^, fi,u = 0, 1, 2, • • • , 2M 



(B6) 

where £^ +) is given by Eq. ([AT]) , g = {gi(a}} are 2M 
arbitrary gauge functions and 



A=(A ,A 1 ,--- ,A 2M ), A =QV, Aj 



A 



(+) 



D = BB T , D = H T : = = f if 



-gtBjk 

(B7) 
(B8) 
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Now, the Ito form [47] of the Langevin equations are 

dQ = ft (vdr + J2 9kdW k ^j , (B9) 



2M 



2M 



da 3 = A$ +) - E 9kB jk )dT + J2 B jk dW k , (BIO) 



fc=i 



fc=i 



where j = 1, 2, • • • , 2M and the Wiener increments dWi 
have the property (dWi(T)dWj(s)) s = 8 13 8(t ~ s)dr 2 [13, 
|49| which can be realized at each dr by real Gaussian 
noises with zero average and variance dr. 



Appendix C: Ito and Stratonovich forms of the 
Langevin equations 



Considering Eqs. (|B9I) and (IB10[) for the Bose- 
Hubbard model, Ito Langevin equations are 

M M M 

dft = ft J ^ UijOLiPj - 1Z^2 n i + Me X! n ' ^ T 



'■3 



i=l 



i=l 



2M 



+ nJ2gkdw k 



fe=i 



da- 7 = 



2i\/ 



Me , 
— rV 

2 



+ i^-aW,, J = 1,2,-.. ,2M. 



(CI) 



r/r 



(C2) 



The Stratonovich stochastic equations are [47J 

2M 

dc^( s ) = (4 s )) dr + Y B^ k dW k: (C3) 
fc=i 

where A^ S) =A fl -^, (j, = 0, 1, 2, • • • , 2M and 

2M 2M 

s p=EE(y-+i;^')i ( c4 ) 

iy = 7=0 



Because B is zero, we can write 



where 



2i\/ 



^=E(Io>+isa)i w - 

2M 2Af 



(C5) 

(C6) 
(C7) 



We have B Q . = Q 9j , a n B Q . = 5j , = 0, 



2 u *7 u ' 



SO 



2M 

So=^Y^ 

3 = 1 



(C8) 



Also 



2 M 2M 



j=i fe=i 

/— 2M 2Af 
2"EE yhjaiPdah - 8kjOe>*d a h'J gj 

j=l k=l 

fu 2M 

k=l 

fu M 

= i^y 2 E + Phdp„) - (a k *d ak * + 0k*d h *)] g, 

fe=i 
fu M 

= ffi V 2 E K a * a °* + " + ^**^»')] 

fe=l 

xn/-K-Ki) = -^E< ( C9 ) 

fc=i 

Considering Eq. (|C8I) . we obtain 

M 

5 = ft^(2. gj 2 - ? ;c/n;') (CIO) 
J'=l 

Furthermore, 9oB = <9o»B = which, according to 
Eq. (IC5|) . gives S^ 1 = 0. Moreover 



2Af 2M 



= -^YYX S ki a j d a k Sjj - Skja 3 * d a k*Sij)a j 
j=i fe=i 

= ~c? (Cll) 

Therefore, we have Si — — ^a 1 . The Stratonovich equa- 
tions, Eq. (|C3|) . are now 

/ 2M \ 

dft (5) = ft Vdr + Y 9kdW k I - y dr 

\ 2M 



ft 



i=i 



" ? 2" n 



fe=i 



271/ 



i=l 



Me , 

2 



(C12) 



S IU 
-fdr + n/ — a J dWi 



J i U . „ ■ 2^ e + £7 
- 2^Wjia - -j(\n>j\ +inj)a J + < 



3 = 1 fe=l 



2 



fir 
(C13) 
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